{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 一、从高斯消元法到 LU 分解\n",
    "\n",
    "在第一节中，学习了矩阵乘法的行视点：初等行变换可以用 $B$ 左乘变换矩阵 $A$ 来表示。\n",
    "在第二节中，讲到线性方程组的求解方法：行简化算法，其实就是对矩阵做一系列的初等行变换。\n",
    "\n",
    "很自然的，行简化算法，能不能用左乘变换矩阵来表示？"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "把 $矩阵A$ 变换成 $上三角矩阵U$ 的方法，又称为高斯消元法。它实际上是对矩阵A做了一系列初等行变换，若用 $变换矩阵E$ 表示这个变换，有\n",
    "\n",
    "$$EA=U$$\n",
    "\n",
    "我们假设有矩阵 $F$ 使 $FE=I$ 成立，方程左乘 $F$，则有\n",
    "\n",
    "$$A=FU$$\n",
    "\n",
    "可以得出的结论是，矩阵 $F$ 是一个下三角矩阵（Lower triangular matrix），一般用 $L$ 表示，于是知道 $矩阵A$可以分解为 $LU$ 的形式。这就叫 LU分解。\n",
    "\n",
    "**用途**：LU分解常用于求解线性方程组，尤其是在我们需要对各种各样的 $b$ 求解方程组 $Ax=b$ 时，LU分解能节省大量的计算步骤。\n",
    "\n",
    "### NOTE\n",
    "\n",
    "1. 注意矩阵A可以为非方阵，若 A 是 $m × n$ 的，则 $L$ 应该是 $m × m$的方阵（变换矩阵E必定为方阵，它的逆L当然也是方阵了），而 $U$ 的形状和 $A$ 相同，为 $ m × n$.\n",
    "1. 从第一条可知，当 $m < n$ 时，方阵 $L$ 会比较大，这不利于计算。但是实际上我们使用的算法得到的 $L$ 是长方形矩阵，而且也完全不需要求 $E$ 的逆。\n",
    "1. 通过上述方法得到的矩阵 $L$ 的对角线元素都为 1."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 举例\n",
    "\n",
    "从矩阵乘法的行视点来看，对矩阵\n",
    "$\\begin{pmatrix}\n",
    "     1 & 3 & 1 \\\\\n",
    "    1 & 1 & -1 \\\\\n",
    "    3 & 11 & 6\n",
    "\\end{pmatrix}$ 做行变换得到 $下三角矩阵U$ 可表示为\n",
    "$$\n",
    "\\begin{pmatrix}\n",
    "     1 & 0 & 0 \\\\\n",
    "     0 & 1 & 0 \\\\\n",
    "     0 & 1 & 1\n",
    "\\end{pmatrix}\n",
    "\\begin{pmatrix}\n",
    "     1 & 0 & 0 \\\\\n",
    "    -1 & 1 & 0 \\\\\n",
    "    -3 & 0 & 1\n",
    "\\end{pmatrix}\n",
    "\\begin{pmatrix}\n",
    "     1 & 3 & 1 \\\\\n",
    "    1 & 1 & -1 \\\\\n",
    "    3 & 11 & 6\n",
    "\\end{pmatrix} = \n",
    "\\begin{pmatrix}\n",
    "     1 & 3 & 1 \\\\\n",
    "     0 & -2 & 2 \\\\\n",
    "     0 & 0 & 1\n",
    "\\end{pmatrix}\n",
    "$$\n",
    "也就是\n",
    "$$\n",
    "E_2 E_1 A = U\n",
    "$$\n",
    "\n",
    "这样就得到了 $U$。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "可以看出 $E$ 也是下三角矩阵： \n",
    "$$\n",
    "E = E_2 E_1 = \n",
    "\\begin{pmatrix}\n",
    "     1 & 0 & 0 \\\\\n",
    "     -1 & 1 & 0 \\\\\n",
    "     -4 & 1 & 1\n",
    "\\end{pmatrix}\n",
    "$$\n",
    "\n",
    "接下来要把 $EA=U$ 变换成 $A=LU$ 的形式，也就是求取使 $LE=I$ 成立的 $L$。\n",
    "显然 $L$ 可以看成是从 $E$ 到 $I$ 的行变换矩阵，应该为\n",
    "$$\n",
    "\\begin{pmatrix}\n",
    "     1 & 0 & 0 \\\\\n",
    "     1 & 1 & 0 \\\\\n",
    "     3 & -1 & 1\n",
    "\\end{pmatrix}\n",
    "$$\n",
    "\n",
    "恰好是个下三角矩阵～"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "所以 $A=LU$ 就是\n",
    "\n",
    "$$\n",
    "\\begin{pmatrix}\n",
    "     1 & 3 & 1 \\\\\n",
    "    1 & 1 & -1 \\\\\n",
    "    3 & 11 & 6\n",
    "\\end{pmatrix} = \n",
    "\\begin{pmatrix}\n",
    "     1 & 0 & 0 \\\\\n",
    "     1 & 1 & 0 \\\\\n",
    "     3 & -1 & 1\n",
    "\\end{pmatrix}\n",
    "\\begin{pmatrix}\n",
    "     1 & 3 & 1 \\\\\n",
    "     0 & -2 & 2 \\\\\n",
    "     0 & 0 & 1\n",
    "\\end{pmatrix}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "下面用 julia 标准库的函数来验证一下"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "LinearAlgebra.LU{Float64,Array{Float64,2}}\n",
       "L factor:\n",
       "3×3 Array{Float64,2}:\n",
       " 1.0   0.0  0.0\n",
       " 1.0   1.0  0.0\n",
       " 3.0  -1.0  1.0\n",
       "U factor:\n",
       "3×3 Array{Float64,2}:\n",
       " 1.0   3.0   1.0\n",
       " 0.0  -2.0  -2.0\n",
       " 0.0   0.0   1.0"
      ]
     },
     "execution_count": 38,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import LinearAlgebra\n",
    "\n",
    "A = [1 3  1\n",
    "     1 1 -1\n",
    "     3 11 6]\n",
    "\n",
    "L, U = LinearAlgebra.lu(A, Val(false))  # LU分解"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "得到 $L$ 和 $U$ 之后，矩阵方程 $Ax=b$ 对任意的 $b$，都有 $LUx=b$，令 $Ux=y$，有\n",
    "$$\n",
    "\\begin{cases}\n",
    "    Ux=y \\\\\n",
    "    Ly=b\n",
    "\\end{cases}\n",
    "$$\n",
    "于是解 $Ax=b$ 就变成了解上述两个矩阵方程。乍一看好像要解的方程还变多了，可仔细想想，$L$ 和 $U$ 都已经是三角矩阵了，直接用回代法就能得到 solution!\n",
    "\n",
    "对任意的 $b$，回代即可得到方程的解，$b$ 越多，节省的计算量就越大！"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 二、矩阵的逆（inverse）\n",
    "\n",
    ">这里**只讨论方阵的逆**，非方阵只可能存在单侧逆。这里不讨论单侧逆，只是简单地认为非方阵不可逆。\n",
    "\n",
    "定义：有 $n × n$ 的 $矩阵A$，若存在 $n × n$ 的 $矩阵C$ 使\n",
    "$$\n",
    "AC = CA = I_n\n",
    "$$\n",
    "\n",
    "成立，则称 $C$ 为 $A$ 的逆，记做 $A^{-1}$。\n",
    "\n",
    ">P.S. 这里只讨论方阵的逆，非方阵的逆很少用到，而且更复杂，基本没见哪本教材谈过。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "不可逆矩阵也称为 **奇异矩阵（singular metrix）**，可逆矩阵称作**非奇异矩阵（non-singular matrix）**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.1 使用矩阵的逆求 $Ax=b$ 的解\n",
    "\n",
    "如果 $矩阵A$ 可逆，有\n",
    "$\n",
    "\\boxed{Ax=b} \n",
    "\\to \\boxed{A^{-1}Ax=A^{-1}b}\n",
    "\\to \\boxed{x = A^{-1}b}\n",
    "$\n",
    "\n",
    "可以看到，$A^{-1}b$ 就是方程的解。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.2 使用矩阵的逆求 LU分解\n",
    "\n",
    "前面我们已经得出了 $$EA=U \\\\ E=E_2 E_1$$\n",
    "\n",
    "于是可以得出 $$A= E^{-1} U = E_1^{-1}E_2^{-1}U$$\n",
    "于是得到 $$L=E_1^{-1}E_2^{-1}$$\n",
    "\n",
    "$E_1$ 和 $E_2$ 都是初等矩阵，都可逆，因此只要有 $EA=U$，那 $L$ 一定存在。\n",
    "\n",
    ">P.S. 这只是理论分析，数值计算中有更好的方法来求 矩阵L。\n",
    "\n",
    "**NOTE**: **永远不要去求矩阵的逆**，矩阵的逆在理论分析上很有价值，但是在数值计算中，效率低下，基本没人用。解矩阵方程一般都用 LU分解。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.2 Gauss–Jordan 求逆法\n",
    "\n",
    "这是一种手算方法，权作了解。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "该方法首先将 $A$ 与 相同大小的单位阵 $I$ 拼接成增广矩阵\n",
    "$\n",
    "\\left(\n",
    "\\begin{array}{c|c}\n",
    "A & I\n",
    "\\end{array}\n",
    "\\right)\n",
    "$，然后对该增广矩阵进行初等行变换，把左边的 $A$ 变换成 $I$，这个时候，右边的矩阵就是 $A^{-1}$\n",
    "\n",
    "该方法的原理也很简单，把 $A$ 变换成 $I$ 可以写成如下形式\n",
    "$$\n",
    "...E_2E_1A=I\n",
    "$$\n",
    "\n",
    "由矩阵的逆的定义 $A^{-1}A=I$ 可知 $A^{-1}=...E_2E_1$\n",
    "\n",
    "在对增广矩阵做初等行变换时，不仅 $A$ 做了该变换，右边的 $I$ 也做了同样的变换，于是右边应该变成了 $...E_2E_1I=A^{-1}I=A^{-1}$，因此右边就是矩阵的逆。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### NOTE\n",
    "\n",
    "1. **可逆的条件**: 从上面可以看出，矩阵A 的逆矩阵 $A^{-1}$ 要存在，就是说 $A^{-1}=...E_2E_1$ 要存在，也就是说**要能通过一系列的初等行变换，把 $A$ 变换成 $I$，即矩阵 $A$ 行等价于单位矩阵 $I$**。**首先把 $A$ 化成阶梯形式，如果它有全为零的行**，这个矩阵就不可能通过对 $I$ 做行变换得到（行变换的数乘要求系数不为0），所以它**一定不可逆**。\n",
    "1. 这里只讨论方阵的逆，方阵化为阶梯形式后，若不存在全为0的行，那它对角线上的元素一定都不为0！也就是说**$n × n$ 的矩阵 $A$ 可逆，等价于它有n个主元列，或者说 $A$ 不存在自由变量。**\n",
    "1. 也可以这样拼接：$\n",
    "\\begin{pmatrix}\n",
    "A \\\\\n",
    "I\n",
    "\\end{pmatrix}\n",
    "$，然后做初等列变换，也能得到同样的 $A^{-1}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 1.0.1",
   "language": "julia",
   "name": "julia-1.0"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.0.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
